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§1.  Introduction 


In  this  paper,  we  present  the  results  of  a  careful  and  detailed  numerical  study  of  the 
small  scale  structures  in  two-dimensional  Boussinesq  convection  in  the  absence  of  viscous 
effects.  In  particular,  we  address  the  issue  of  whether  a  finite  time  singularity  can  form  from 
a  smooth  initial  data.  Recently  in  their  numerical  study  of  the  same  problemt^'*^,  Pumir  and 
Siggia  observed  that  the  cap  of  a  symmetric  rising  bubble  (with  smooth  density  variation) 
collapses  in  a  finite  time.  In  contrast  with  their  observations,  our  results  suggest  that  such 
collapse  cannot  occur  in  a  finite  time.  The  strain  rate  associated  with  the  intensification 
of  the  density  gradient  saturates,  implying  an  exponential  decay  for  the  thickness  of  the 
iront.  This  is  reminiscent  of  the  situation  found  in  vortex  reconnectiont'’'^’*^^;  when  two 
vortex  tubes  are  brought  together,  the  axial  strain  rate  saturates  and  the  core  of  the  tubes 
undergoes  enormous  deformation  to  avoid  reconnection  (in  the  absence  of  viscosity).  Thus 
the  inviscid  solution  manages  to  escape  from  forming  a  finite  time  singularity. 

There  are  two  main  motivations  for  the  study  of  the  two-dimensional  Boussinesq  con¬ 
vection.  One  comes  from  the  jjotential  relevance  of  this  problem  to  the  study  of  atmospheric 
and  oceanographic  turbulence,  as  well  as  other  astrophysical  situations  where  rotation  and 
stratification  play  a  dominant  role.  The  second  motivation  comes  from  the  fact  that  from 
a  computational  viewpoint,  this  is  the  simplest  among  the  class  of  incompressible  flows 
which  exhibit  vorticity  intensification.  In  particular,  it  is  an  open  question  whether  the 
l)aroclinic  generation  of  vorticity  leads  to  a  finite  time  singularity.  Such  singularities,  if  ex¬ 
ist,  provide  an  effective  mechanism  for  the  cascade  of  energy  to  small  scales.  This  scenario 
also  provides  a  convenient  basis  for  various  turbulence  theories  which  assume,  in  one  form 
or  another,  that  the  (ensemble  average  of)  rate  of  viscous  dissipation  of  energy  remains 
finite  in  the  limit  of  vanishing  viscosity,  implying  the  occurrence  of  finite  time  singularities 
in  many  Euler  flowst’^l. 

There  is  a  well-known  analogy  between  the  two-dimensional  Boussinesq  convection 
and  the  three-dimensional  axisymmetric  swirling  flows: 

f  Vt  +  tlVr  +  IUV2  -t-  ytxw  =  0 

'  '  '  (  -b  ULVr  +  ^  UU!  —  j:(v^)z  =  0 
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Here  u  =  tx  er  +  e  e^+u>  is  the  velocity,  u}  =  u^—Wr  is  the  azimnthal  vorticity.  Comparing 
(1.1)  with  (2.1)  -  (2.4),  we  see  that  the  centrifugal  force  plays  a  role  similar  to  gravity,  and 
the  azimuthal  circulation  plays  a  role  similar  to  density.  Grauer  and  SiderisI®]  were  the 
first  to  seek  finite  time  singularities  in  this  restricted  class.  Although  their  numerical  result 
still  remains  inconclusive,  it  has  stinmlated  a  lot  of  recent  work  along  the  same  direction, 
including  the  present  study. 

A  related  problem  was  studied  earlier  by  Childress^^^  where  he  considered  nearly  two- 
dimensional  Euler  flows  and  derived  effective  equations  governing  their  slow  variation  using 
contour  averaging  methods.  Under  a  change  of  variables,  the  effective  equation  takes  a  form 
similar  to  the  axisymmetric  Euler  equations  with  nonstandard  connection  between  circula¬ 
tion,  radius  and  the  angular  velocity  component.  Childress  went  further  to  study  numeri¬ 
cally  a  simplified  version  of  his  effective  equations  and  observed  finite  time  singularities  for 
the  simplified  model.  This  was  interpreted  as  a  signal  for  the  re-three-dimensionalization 
of  the  original  (nearly  two-dimensional)  flow. 

This  paper  is  organized  as  follows.  In  §2.  we  formulate  the  problem  and  present 
some  preliminary  mathematical  remarks.  In  §3.  we  describe  the  numerical  methods  used 
to  study  this  problem.  The  numerical  results  are  presented  in  §4.  Some  discussions  and 
concluding  remarks  are  made  in  §5. 


§2.  Formulation  of  the  Problem; 


The  equations  describing  Boussinesq  convection  are  the  following: 

p/  -h  u  ■  Vp  =  0 

(2.1)  ^  u,  u  •  Vu  -h  Vp  =  -  (°) 

V  ■  u  =  0 


Here  p  is  the  density  (this  should  be  the  temperature  and  denoted  by  9  ov  T.  But  we  are 
accustomed  to  calling  it  density,  and  therefore  denote  it  by  p)  ,  u  =  (w,e)  is  the  velocity, 
p  is  the  pressure.  We  have  normalized  the  gravitational  constant  to  be  1. 
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It  is  convenient  to  write  (2.1)  in  the  streamfunction  -  vorticity  formulation: 


r  +  u  •  V/j  =  0 

(2.2)  <  tJt  +  u  •  Vcj  = 

I  —  All’  =  uj 

Here  io  is  the  vorticity,  xj)  is  the  stream  function: 

(2.3)  UJ  —  Vx  -  Uy  u  =  il>y,  V  =  —il’x- 
Introducing  the  material  derivative  Dt  =  dt  +  u.V,  (2.2)  becomes  simply 

(2.4)  Dtp  =  0,  Dtuj  =  -  px. 


It  is  straightforward  to  show  that  if  the  initial  data  u(a;,y,0)  =  Uo{x,y)  and  p{x,y,0) 
=  po{x,y)  are  smooth  enough,  then  the  solution  to  (2.1)  exi.«ts  and  remains  smooth  for  a 
short  time.  It  is  not  clear  whether  any  solution  would  loss  its  regularity  at  a  finite  time. 
However,  following  Beale,  Kato  and  Majdat^l,  one  can  show  that  if  a  solution  losses  its 
regularity,  then  the  density  gradient  has  to  blow  up.  More  precisely,  we  have: 


Theorem:  Define  the  norms: 


( 


ll/(- 


S  / 

•  t»l  ,«2^0  ^2 

\  ai  +a2<m 


^ai+rt2  j 


\ 


1/2 


dx°‘’ 


dx  dy 


[/(Oloo  =  ^  max  \f{x,y)\. 

Assume  that  for  some  m  >  2,  ||u(-,0)||m  +  ll/^(’5  0)||m  is  finite,  but  there  exists  a  T*  such 
that  ||u(-,T*)||m  +  l|p(-,T*)|lm  =  +00.  Then 

rT' 


(2.5) 
and 

(2.6) 


/ 

Jo 


t)\oodt  =  +00 


/  /  \pxi-,s)\oo  dsdt  —  +00. 

'0  Jo 


Proof:  We  will  only  give  a  sketch  of  the  proof  since  it  follows  closely  the  argument 
of  Beale,  Kato  and  Majdat^l.  We  will  use  C  to  denote  a  generic  constant. 
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Straightforward  energy  estimates  lead  to 

<  |Vu(.,<)U  U-Ml, 

<  lVu(-,()U  ||u(-,t)£  +  l|p(-,t)£. 

The  logarithmic  Sobolev  inequality  gives  us,  for  m  >  2 

|Vu(-,t)U  <  CK-,t)|oo(l  +  log  ||u(.,<)l|m). 


Therefore,  if  we  let 

y{t)  =  ||/)(-,<)!lm  + 

we  have 

y{t)  <  C|u;(-,<)joo  (1  +  log  y{t))y{t). 
Solving  this  differential  inequality,  we  get 

y(t)  <  y(0). 

On  the  other  hand,  since  we  have 


l^('>^)loo  ^  lP*(')f)loo- 


Hence 


and 


|c<;(-,  t)l<x)  <  f  |pi(', 5)100  ds  C 

Jo 


c  f‘  f’  lPx(  .liloodedf 

V{t)  <  y(0). 

Since  y(0)  is  finite,  we  conclude  that  if  y(T*)  =  +00,  then  (2.5)  and  (2.6)  hold. 


§3.  The  Numerical  Methods; 

Computing  singular  or  nearly  singular  solutions  is  a  difficult  task,  particularly  so  in 
the  context  of  incompressible  flows.  To  obtain  maximum  information  one  has  to  push  the 
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mimerical  method  to  the  point  where  the  flow  is  only  marginally  resolved.  In  this  situation 
any  numerical  method  is  likely  to  exhibit  its  own  artifact.  Such  numerical  artifacts  may 
not  go  away  under  simple  mesh  refinement  checks,  and  they  do  not  necessarily  manifest 
themselves  in  scales  comparable  to  the  grid  size.  Therefore,  there  is  a  real  danger  of  being 
mislead  by  a  particular  numerical  result.  To  avoid  this,  we  have  solved  (2.1)  using  two 
different  numerical  methods:  a  Fourier  -  collocation  method  and  an  ENO  scheme.  When 
the  solution  develops  large  gradients,  the  Fourier  -  collocation  method  usually  exaggerates 
the  situation,  whereas  the  ENO  scheme  tends  to  smear  out  the  large  gradients.  By  com¬ 
paring  the  numerical  results  using  the  two  different  methods,  we  can  Ijetter  judge  what 
phenomena  is  likely  to  be  physical.  Below  we  will  describe  separately  the  two  methods. 

3.1  Spectral  Method  with  Smoothing 

This  is  the  standard  Fourier-collocation  methodt®^  with  smoothing  or  de-aliasing. 
Roughly  speaking,  the  differentiation  operator  is  approximated  in  the  Fourier  space,  while 
the  nonlinear  operations  such  as  multiplications  are  done  in  the  physical  space.  We  used 
the  intrinsic  Cray  FFT  routines  which  considerably  enhanced  the  performance  of  the  code. 

Since  the  solutions  have  large  gradients,  it  is  crucial  to  add  filters  to  the  spectral 
method  in  order  that  the  numerical  solutions  do  not  degrade  catastrophically  if  some  part 
of  the  steep  gradients  are  not  adequately  resolved.  A  roljust  way  of  adding  the  filterst'^^l 
is  to  replace  the  Fourier  multiplier  ikj  for  the  differentiation  operator  by 
where 

k 

(3.1)  <^[k)  =  e  ^ -N  ^  ,  for  |A'|  <  iV. 

Here  N  is  the  numerical  cutoff  for  the  Fourier  inodes,  nif  is  the  order  of  the  filter,  and  o 
is  chosen  so  that  ^{N)  =  =  machine  accuracy.  The  machine  accuracy  on  Cray-YMP 

with  single  precision  is  roughly  10“'  b  Denote  by  F  and  respectively  the  forward  and 
backward  Fourier  transform  operators,  then  the  numerical  derivative  is  ev'ahiated  as 

(3.2)  D^f  =  F-\ik^{\k\))Ff. 
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The  accuracy  of  such  an  approximation  scheme  depends  on  the  parameter  m/.  For  smooth 
functions  /(x),  we  have 


(3.3)  ||/'(x)  -  Z)yv/(x)||  =  0(iV-'"0 

We  will  denote  the  Fourier-collocation  method  with  m/-th  order  filter  as  SPmf.  Unless 
otherwise  stated,  the  results  presented  in  §  4.2  were  computed  with  rrif  =  10. 

3.2  The  ENO  Scheme 

The  ENO  (Essentially  Non-Oscillatory)  scheme  is  used  for  the  convection  (spatial)  part 
of  the  flow.  We  use  ENO  schemes  based  on  point  values  and  numerical  fluxes,  developed 
by  Shu  and  Osher  in  [16],  [17];  see  also  [18]. 

To  apply  ENO  approximations,  the  equation  is  first  written  in  a  conservation  form. 
For  example  the  first  equation  in  (2.2)  is  written  as 


Pt  +  {up)x  +  {vp)y  =  0 


(3.4) 


The  ENO  operator  is  then  applied  to  each  of  the  conservative  derivatives  in  a  dimension- 
by-dimension  fashion;  when  approximating  {up)x,  y  is  fixed.  Unlike  the  compressible  flow, 
the  incompressible  flow  equations  are  naturally  written  in  characteristic  form.  Thus  no 
expensive  characteristic  decomposition  is  needed.  Upwinding  can  be  simply  determined 
by  the  signs  of  u  and  v. 


We  shall  only  briefly  describe  the  approximation  of  a  .single  derivative,  say  More 
details  can  be  found  in  [16],  [17],  and  [18]. 


The  conservative  approximation  is  of  the  form: 

(/.)>==  (3.5) 

which,  for  r-th  order  ENO  scheme,  approximates  the  derivative  to  r-th  order: 

i  +0(Ax')  (3.6) 
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where  the  numerical  flnx  fj^  ^  is  obtained  by  interpolating  the  point  values  of  /  on  a  stencil 
of  ?■  +  1  consecutive  grid  points.  The  stencil  is  chosen  inductively  as  follows.  For  ?•  =  1, 
we  choose  the  stencil  to  be  [j-1,  j]  or  [j,  j+l]  depending  on  the  sign  of  ii.  For  r  >  1,  left  or 
right  neighboring  points  are  added  to  the  stencil  at  the  previous  level  r  —  1  according  to 
the  absolute  value  of  the  divided  differences  they  each  give.  In  most  of  our  calculations, 
we  used  r  =  3.  We  will  denote  this  method  by  ENOS.  Notice  that  the  scheme  is  actually 
1'  T  Tth  order  in  L\  norm.  That  is,  the  third  order  ENO  scheme  we  use  is  actually  fourth 
order  accurate  in  the  L\  norm. 

The  potential  equation  in  (2.2)  is  solved  with  a  fourth  order  central  differencing  (im¬ 
plemented  in  the  Fourier  space  via  FFT)  plus  a  fourth  order  exponential  filtering  in  the 
Fourier  space  described  above.  This  guarantees  fourth  order  accuracy  and  in  most  cases 
avoids  instability.  We  are  not  sure  whether  ENO  should  or  can  be  applied  to  this  potential 
solver,  since  it  is  an  elliptic  equation  with  possibly  a  singular  right-hand-side.  Unlike  the 
compres.sible  case,  some  small  oscillations  can  still  be  seen  in  the  ENO  solution,  probably 
due  to  this  potential  solver. 

For  the  temporal  discretization,  we  u.sed  Runge-Kutta  methods  of  various  order  de¬ 
signed  inl'^1.  No  major  difference  between  the  3rd,  4th  and  5th  order  methods  were  found 
in  the  numerical  results.  It  seems  to  be  a  general  fact  that  temporal  accuracy  is  much  less 
important  than  the  spatial  accuracy.  We  u.sed  the  third  order  version  most  often  since  it 
only  recpiires  three  auxiliary  arrays,  whereas  the  fourth  order  version  requires  five  auxiliary 
arrays.  We  take  initial  data  that  is  periodic  with  period  D  where  D  =  [0,27r]  x  [0,27r].  The 
results  reported  below  were  computed  iising  CFL  equal  to  0.5.  This  is  very  much  within 
the  stability  region  of  these  methods. 

We  point  out  that  our  experience  strongly  favors  the  use  of  high  order  schemes  l)e- 
cause  of  their  small  dispersive  errors.  In  a  marginally  resolved  situation,  even  though  the 
filtered  spectral  method  does  generate  small  oscillations  near  the  front  wdiere  the  density 
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experiences  a  jump,  the  oscillations  are  effectively  localized  near  the  front.  This  is  l)ecaiise 
the  spurious  numerical  oscillations  travel  at  roughly  the  right  speed  which  is  the  sjjeed  of 
the  front.  The  higlu'r  the  order  of  the  filter,  the  more  accurate  the  proi)agation  sj)eed.  the 
more  localizetl  the  oscillations.  Of  course  the  Gibbs  phenomena  eventually  ])revails  and 
the  numerical  result  beconu's  m^ise. 


§4.  Numerical  Results 

4.1  Formation  of  the  front 


We  have  numerically  integrated  (2.1)  with  a  variety  of  initial  data.  The  initial  density 
is  usually  a  j)erturbation  of  the  uniform  state,  and  the  perturbations  are  localized  in  each 
period.  The  early  time  evolution  of  s\ich  flows  is  characterized  by  the  formation  of  a  front 
ac'ross  which  density  varies  sharply.  As  time  evolves,  the  front  gets  incieasingly  sharper 
and  th('  tail  starts  to  roll  up.  Figures  1-3  present  the  time  development  of  such  events 
for  the  initial  data 


(4.1) 


cj(r,y,0)  =  0 

p(j-,,(/,0)  =  50pi(x,y)p2(j',y)(l  -  di(-r.y)) 


where 


I  exp  (l  -  if.r^ +  (j/-7r)2  <  TT^ 

[  0  otherwise. 


p-A-f 


y)  r=  (  ^-xp  (l  -  ,  if  1-^  -  27r|  <  I.DStt 


I  0  otherwise 

Figure  1  is  the  density  contour  at  f  =  1.6  computed  using  SPlO  on  a  512^  grith  At  this 

time,  the  flow  looks  roughly  like  a  rising  b\ibble.  Figure  2  presents  the  same  information 
at  /  =  2.5.  By  now  the  outer  boundary  of  the  bubble  has  become  a  sharp  front.  We  notice 
that  as  the  bubble  rises,  it  leaves  behind  a  long  and  thin  filament  of  light  fluid.  This  is  a 
check  on  the  amount  of  numerical  diffusion  present  in  the  scheme.  A  low  order  method 
with  numerical  viscosity  (needed  to  stabilize  the  front)  will  destroy  the  thin  filament.  In 
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Fillin'  3a  and  3'*  we  })ieseiit  respect iv('ly.  the  coiitoiirs  of  density  and  vorticity  at  /  =  2.9. 
Now  tilt'  filanu'iit  of  light  Hnirl  is  thinner  and  longer,  and  the  trailing  edg(>  of  the  front 
has  rolh'd  nj).  From  Figiire  3b  we  sc'e  that  most  of  the  vorticity  resides  tm  the  front,  the 
rt'st  of  the  flow  domain  has  very  little  vorticity.  This  fact  is  more  drastically  displayed  in 
Figure  4a  and  41)  when'  we  plot  resjiectively.  the  slice  of  u,'  and  pj.  at  y  =  Since  a,'  and 
~l>j.  have  the  same  sign,  from  the  vorticity  ecpiation  (2.2)  we  see  that  the  vorticity  peaks 
will  he  incn-asingly  and  monotonically  .sharper.  However,  it  a])])ears  that  this  inechai  ism 
of  l)aroclinic  vorticity  generation  only  giv('s  rise  to  an  ex])onential  growth  of  vorticity.  We 
will  come  hack  to  this  point  later. 

To  better  a])pr('ciate  the  nearly  singular  behavior  and  the  computational  difficulties 
W('  fact'.  W('  will  show  a  ft'w  more  plots  for  the  profih's  of  the  density  and  the  velocity. 
Figurt'  5a  displays  the  evolution  of  the  density  along  the  symmt'try  axis  ./■  =  ~  at  times 
f  =  0.5.  1.  1.5.  2  and  2.5.  It  shows  clearly  the  formation  of  a  front,  similar  to  the  formation 
of  shock  fronts  in  the  .solutions  of  the  Burgers  etjuafion.  .4t  t  =  2.5.  SPIO  on  the  512“^  grid  is 
not  rt'solving  tlu'  flow  and  small  nuim'rical  oscillations  appear  on  the  profile.  Howevc'r  the 
numerical  oscillations  die  out  on  the  n'fined  grid  1024'.  The  result  of  the  latter  calculation 
is  sup«'rimposed  in  Figure  5a  on  th<'  solution  of  the  512'^  Sri<^l-  hi  Figure  5b  and  5(  we 
show  the  slices  of  p  and  c  respectively,  at  y  =  tt.  f  =  2.5.  From  these  graphs  one  might  l)e 
tempt<’fl  to  conclude  that  jumps  hav('  occurred  in  p  and  cusps  have  occurred  in  v.  However, 
a  clos(’r  look  suggests  that  th('  cuspdike  behavior  in  r  and  the  discontinuotis  behavior  in 
p  go  away  under  refinement.  In  Figure  5d  we  coinpan'  the  numerical  results  for  r  on  the 
512"^  and  1024^  grids  in  the  cusp-like  regions.  It  is  seen  that  on  the  1024“^  grid  the  solution 
looks  much  less  singidar.  Similar  phenomena  were  observed  for  p. 

These  restdts.  particularly  Figure  2  an<l  3a  motivate  the  following  question;  Is  ' 
possible  that  the  solution  losses  its  regularity  at  the  satire  time  at  the  e'litire  front?  For 
this  to  happen  there  are  only  two  possibilities.  One  is  that  while  r  develojrs  cusjrs.  p  stays 
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smooth  across  the  front.  This  is  rnletl  out  by  Theorem  1  since  it  dictates  that  p  and  uj 
must  become  singidar  at  the  same  time.  The  second  possibility  is  that  p  develops  jumps 
across  the  emtire  front.  In  this  case  it  is  necessary  that  fluids  with  intermediate  densities 
are  suddenly  swept  to  the  back  of  the  bubble.  Clearly  this  requires  infinite  I'elocity  whereas 
an  of  our  numerical  restdts  give  velocities  that  are  very  well  bounded. 

Before  ending  this  subsection,  we  present  the  comparison  of  the  numerical  results 
computed  using  the  two  different  methods  describetl  in  §3.  Figure  6a  and  Gb  display 
re.-,j)ectively,  the  time  history  of  the  maximum  and  minimum  density  computed  using  SPIO 
on  256'^.  612^,  and  1024^  grids  and  EN03  on  the  512“^  grid.  These  quantities  should  be 
conserved  by  the  exact  solutions.  As  expected,  EISO  does  a  much  better  job  in  avoiding 
overshoots  and  undershoots.  Figure  6c  compares  the  numerical  results  for  p  at  ?/  =  tt,  t  = 
2.5,  computed  using  SPIO  and  EN03  on  512^  grid.  The  result  of  SPIO  undershoots  much 
more  than  the  result  of  EN03,  although  the  latter  also  contains  some  small  numerical 
oscillations.  On  the  other  hand,  although  not  shown  here,  ENO  is  usually  less  accurate 
than  the  Fourier-collocation  method  and  contains  more  numerical  dissipation. 

4.2  Evolutions  of  the  bubble  cap  and  other  aspects  of  the  flow. 

When  the  bubble  rises,  lighter  fluid  has  larger  acceleration  and  heavier  fluid  has  smaller 
acceleration.  This  results  in  the  formation  of  fronts.  However,  once  the  front  is  formed, 
the  pressure  gradient  yccomes  iniijortant  across  the  front.  In  fact,  as  is  shown  in  Figures 
7a  and  7b.  the  formation  of  increasingly  larger  pressure  gradients  may  reverse  the  initial 
picture.  In  Figures  7a  and  7b.  we  plot  the  slice  of  —py.  p,  and  the  acceleration  — (p  +  Py) 
along  the  symmetry  axis  =  tt  at  two  different  times  i  =  0.5  and  f  =  2.  At  f  =  0.5,  the 
acceleration  —{p^  Py )  decreases  acro.ss  the  front,  whereas  at  f  =  2  it  increases  across  the 
front.  (The  density  p  always  decreases  acro.ss  the  front,  see  Figure  5a)  This  implies  that 
the  velocity  difference  bc'tween  the  fluitl  particles  at  the  tip  and  the  back  of  the  front  will 
not  increase  further  after  f  =  2,  assuming  that  the  picture  remains  valid. 
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From  another  viewpoint,  the  fluid  on  the  symmetry  axis  satisfies 

(  pt  +  vpy  =  Q 

)  j  =  0 

I  Vt  +  VVy  =  -  (p  +  Py). 


Let  ^  =  py^  7]  =  Vy.  Then  we  have 


6  +  V^y  ^  -Tj^. 


The  growth  of  ^  is  controlled  by  the  strain  rate  t].  In  Figure  8a  we  exhibit  the  time  history 
of  the  maximum  of  |7/|.  During  the  time  interval  displayed,  the  maximum  of  |7;|  always 
occurs  at  the  bubble  cap.  The  results  of  several  different  calculations  agree  remarkably 
well;  the  strain  rate  7/  saturates  at  about  t  =  1.7.  Con.sequently,  the  maximum  of  |^j  can 
at  most  grow  exponentially.  This  is  indeed  so,  as  can  be  seen  from  Figure  8b.  Other 
quantities,  such  as  w  and  pi,  also  grow  at  an  exponential  rate. 

We  have  also  done  calculations  with  other  .sets  of  data.  Here  we  will  briefly  report  the 
results  of  a  few  cases. 

Figures  9a  -  9c  display  the  results  for  initial  data 

u;(x,j/,0)  =  0 


p{x,y,0)  =  p'iix.y)  = 


(277  — j:)e  2 


if(x  —  tt)^  +  (y  —  <  2.5^ 

otherwise. 


The  contour  curves  of  p  at  t  =  4.5  and  t  =  5  are  shown  respectively  in  Figures  9a  and 
9b,  whereas  Figure  9c  shows  the  time  history  of  the  strain  rate  [t;!  at  the  bubble  cap.  We 
observe  that  the  strain  rate  begins  to  saturate  at  t  ^  4.25  and  then  increases  quickly  at 
t  «  4.8.  This  corresponds  to  the  time  when  a  mushroom  forms  out  of  the  bubble  cap,  and 
the  front  changes  its  global  configuration.  We  see  that  eventually  the  strain  rate  saturates 
again  at  about  t  =  5.25.  We  expect  that  the  strain  rate  at  the  cap  will  not  increase  unless 
the  flow  changes  its  geometry  at  the  cap. 
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Figiin's  10a  10c  present  the  numerical  results  for  a  set  of  more  complicated  initial 

data.  The  flata  is  designed  so  that  at  lat<'r  times  the  bubble  front  will  develop  two  bumjis 
instead  of  one  as  in  Figure  9b.  Again  the  increa.se  of  the  strain  rate  at  the  cap  is  accom¬ 
panied  by  a  change  of  the'  front  configiiration:  after  t  —  3.5.  the  front  is  no  longe'r  conve'x 
(as  it  was  liefore  /  =  3.5).  We  were'  not  able  to  carry  this  computation  furthe'r  to  see  the' 
eve'iitual  saturatie)n,  since  the  siele  arms  of  the  bubble  become  so  thin  that  tlu'y  cannot  be 
re'solved  on  a  1024"  grid.  It  is  not  cle'ar  whetlu'r  a  pinch  actually  occurs  at  the'  siele  arms. 


§5.  Discussions  and  Conclusions 

Our  numerie-al  re'sults  strongly  sugge'st  that  the  collapse  e)f  the  bulible  cap  observed  l)y 
Pumir  and  Siggiaf’’*^  is  unlikely  te)  ea-cur  in  a  resolved  calculation.  The  increase  of  the  strain 
rate  is  always  associate'd  with  a  global  change  of  the  flow  character.  The  contribution  of  the 
local  strain  tends  to  saturate,  The  global  character  of  the  flow  is  impen'tant.  Conseque’iitly 
any  numerical  proceelure'  which  e'lnploys  se)me  sort  of  local  approximation  is  open  to  serious 
suspicion.  This  inchtdt's  the  ones  presented  in^'^K  However,  we  did  not  address  the  issue 
wlu'ther  other  tyjtes  of  singularities  can  occur. 

Let  us  reiterate  the  argument  presented  in  §4.2.  Along  the  symmetry  axis  .r  =  we 

have 


(5.1) 


(  D,^  =  -  tj(, 

{Dili  =  -{()y  +  Pyy)  -  l}' 


wIk'H'  (f  =  py.  tj  =  Vy  are  both  negative  at  the  front.  From  (2.1)  we  also  have 


o.Z 


-py  -  Ap  =  Ul  +  I'l  +  2lly  Vr. 


On  the  .symnu'try  axis,  this  reduces  to 


(5.3) 


Py  Pyy  ~  ■ 
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Therefore,  we  have 


(5.4) 


f 

[DtJ]  =  pr 


+  7/2. 


As  long  as  the  front  remains  relatively  smooth,  Pu-  is  negligible  compared  to  if .  Conse¬ 
quently  ^  and  7/  will  stay  bounded.  In  fact  since  near  the  cap  py  <  0,  pyy  <  0,  we  obtain 
a  necessary  condition  for  collapse  at  the  cap: 


(5.5) 


-27/2  ^ 


In  particular,  the  front  itself  has  to  develop  infinite  curvature.  This  is  consistent  with  the 
argument  presented  in  section  5  off'‘*K 

The  key  difference  between  the  computations  reported  int'‘^l  and  the  computations 
reported  here,  is  that  int'"*]  a  high  frequency  instability  (with  a  scale  that  is  comparable  to 
the  thickness  of  the  front)  is  observed  at  the  cap,  whereas  in  all  our  resolved  calcidations 
we  have  not  seen  the  manifestation  of  such  an  instability.  Since  the  front  is  extremely  thin, 
near  the  cap  the  situation  resembles  that  of  the  Rayleigh- Taylor  instability.  Therefore  a 
linear  stability  theory  will  indeed  predict  that  the  front  is  unstable  to  short  wavelength 
perturbations.  However,  the  relevant  question  here  is  not  whether  perturbations  can  grow, 
but  whether  short  wavelength  perturbations  can  exist  in  a  dynamic  context.  In  other 
words,  the  issue  is  whether  the  exact  solutions  of  (2.1)  can  supply  perturbations  with 
wavelength  on  the  scale  of  the  thickness  of  the  front. 

We  want  to  emphasize  that  the  question  we  are  asking  is  a  rather  academic  one.  In 
actual  physical  experiments,  various  external  noises  provide  ample  sources  of  such  pertur¬ 
bations.  The  solution  we  are  .seeking  numerically  may  hardly  occur  in  experiments.  (This 
is  yet  another  argument  for  the  importance  of  scientific  computing:  It  has  the  potential  to 
carry  out  much  more  controlled  experiments.) 

Before  discussing  the  question  raised  above,  let  us  recall  two  problems  that  bear 
some  similarity  with  the  case  we  are  studying.  The  first  is  the  classical  Kelvin-Helmholtz 
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instability  of  a  vortex  sheet.  The  braid  that  connects  two  consecutive  rolls  in  a  Kelvin- 
Hehnholtz  roll-up  is  still  subject  to  Kelvin-Hehnholtz  instability,  although  the  local  vortex 
sheet  strength  is  decreasing  in  time  and  the  braid  is  under  straining.  This  kind  of  instability 
manifests  itself  in  numerical  computations  in  such  a  way  that  the  round-off  error  causes 
small  si)urious  roll-ups  along  the  braids.  This  phenomenon  is  studied  very  carefully  l)y 
Krasnyt’’^  and  is  cured  by  introducing  a  nonlinear  filter  that  keeps  the  numerical  solution 
effectively  analytic. 

The  second  problem  is  the  liehavior  of  an  expanding  liubble  in  a  Hele-Shaw  cell  with 
very  small  surface  tension.  The  presence  of  the  surface  tension  makes  the  interface  between 
two  viscous  fluids  linearly  stable.  Nevertheless,  the  interface  is  nonlinearly  unstable  and 
prefers  to  develop  the  characteristic  viscous  fingers  and  the  tip-splitting  of  these  fingers 
observed  in  experiments.  For  example  the  amount  of  perturbation  needed  to  upset  the 
linear  stability  of  a  viscous  finger  and  cause  it  to  split  is  exponentially  small  in  terms  of 
the  surface  tensions'll  .  In  the  numerical  studies  of  Dai  and  Shelley!®!  ^  the  lower  precision 
calculations  give  rise  to  finger-like  behavior.  One  might  argue  that  such  fingers  are  physical, 
on  the  basis  that  they  are  observed  in  experiments.  But  they  are  not  faithful  solutions 
of  the  original  dynamic  problem.  Indeed  as  was  shown  by  Dai  and  Shelley,  the  initial 
instability  go  away  in  the  high  precision  calculations  (128  bits). 

Retxirning  to  the  problem  we  are  studying,  the  linear  stability  of  the  cap  of  a  rising 
bubble  is  investigated  by  Batchelor  et.  al.!^’*'*!.  In  particular,  Pumir  and  Siggia  show  that 
perturbations  with  wavelength  6  (the  thickness  of  the  front)  can  ami)lify  by  a  factor  of  ^ 
(where  70  is  the  radius  of  curvature  at  the  cap)  if 


(5.6) 


7r 

6 


26s'^ 


2((— 

<  e  2^ 


1) 


Heie  .s  is  the  strain  rate  at  the  cap.  We  remark  that  an  infinitesimal  disturbance  eventually 
decays  since  the  local  strain  increases  its  wavelength  exponentially. 
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Although  in  some  cases  we  have  carried  out  the  computation  well  beyond  the  limit 
set  by  (5.6),  we  have  not  observed  the  kind  of  instabilitj'  at  the  cap  reported  I)y  Pumir 
and  Siggial^^l.  One  such  instance  is  displayed  in  Figure  11  where  we  have  roughly  6  = 

')c  =  =  2.  (5.6)  becomes  2048  <  1.109  x  10“^  and  is  very  well  sati.sfied.  Still  no 

instability  on  the  front  is  observed. 

A  similar  situation  occurs  in  the  classical  Rayleigh- Taylor  problem.  The  interface 
between  a  heavy  fluid  on  top  and  a  light  fluid  at  the  bottom  is  linearly  unstable  and 
typically  evolves  to  mushroom-like  or  plume-like  structures.  The  tips  and  tails  of  the 
structure  are  still  subject  to  linear  instability.  But  no  instability  at  the  tips  or  tails  has 
been  reported,  even  though  the  problem  has  been  extensively  studiedt^'’’^®^. 

We  remark  that  for  the  collapse  to  occur,  a  cascade  of  instabilities  has  to  develop  at 
the  cap  with  increasingly  smaller  temporal  and  spatial  scales.  We  have  argued  that  such 
an  instability  is  unlikely  to  arise  for  the  true  exact  solutions  of  (2.1).  However,  since  the 
front  is  very  sensitive  to  numerical  errors,  any  uncontrolled  truncation  can  easily  trigger 
the  instability  and  lead  to  erroneous  results. 

With  respect  to  the  issue  of  formation  of  small  scale  structures,  our  calculation  sug¬ 
gests  that  the  cascade  of  energy  to  small  scales  is  dominated  by  the  formation  of  thin 
fronts.  Roll-ups  that  are  characteristic  of  two-dimensional  Kelvin-Helmholtz  instability 
also  occur  but  play  a  less  important  role.  Occasionally  we  also  observe  small  rolls  formed 
out  of  big  rolls,  evidence  of  what  is  suggested  in  Richardson's  well-known  quote.  But  these 
rolls  carry  a  very  small  amount  of  energy.  The  scenario  of  a  cascade  through  a  sequence 
of  roll-ups  on  increasingly  smaller  scales,  if  at  all  possible,  has  negligible  contribution  to 
the  overall  cascade  mechanism. 

The  evolution  of  the  energy  spectrum  for  the  first  set  of  data  is  plotted  in  Figure 
12.  There  is  an  apparent  similarity  to  the  spectrum  observed  by  Brachet  et.  al.  for  the 
Taylor-Green  vortex^^i.  The  near  singular  behavior  creates  an  algebraically-decaying  part 
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in  the  biilk  of  the  spectrum.  However  the  spectrum  does  not  seem  to  settle  down  to  any 
stationary  state. 

To  conclude,  let  us  emphasize  that  the  ultimate  answer  to  the  question  addressed  here 
must  be  analytical.  Numerical  studies  can  at  best  provide  partial  evidence.  For  incom¬ 
pressible  flows,  it  has  proven  very  difficult  to  obtain  analytical  results  on  the  regularity  of 
the  solutions.  As  a  result,  considerable  efforts  were  devoted  to  the  numerical  aijproach. 
However,  computing  singular  or  nearly  singular  incompressible  flows  is  also  an  extremely 
difficult  task.  Not  surprisingly,  the  subject  of  seeking  singularities  numerically  for  incom¬ 
pressible  flows  is  full  of  controversy.  The  contradicting  conclusions  reached  in  the  present 
paper  and  the  work  of  Pumir  and  Siggiat'"*^  certainly  adds  one  more  to  the  list  of  many 
unsettled  issues.  It  is  our  intention  that  these  numerical  work  will  generate  sufficient  inter¬ 
est  for  further  analytical  study  and  the  evidence  presented  will  contribute  to  the  ultimate 
resolution  of  the  issue  raised. 
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Figure  1.  Density  contour  at  t=1.6  for  initial  data  (4.1) 
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Figure  2.  Density  contour  at  t=2.5  for  initial  data  (4.1). 


20 


Figure  3a.  Density  contour  at  t=2.9  for  initial  data  (4.1). 
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Figure  3b.  Vorticity  contour  at  t=2.9  for  initial  data  (4. 
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Figure  5a.  Evolution  of  density  along  the  symmetry  axis  x  =  rr  at  f  =  0.5, 1. 1.5. 2. 
and  2.5.  computed  using  SPIO  on  a  512^  grid.  The  result  at  f  =  2.5  computed  on  a  1024 
grid  is  also  superimposed.  Notice  that  the  small  numerical  oscillations  disappear  on  the 
refined  grid. 


Figure  5b.  Slice  of  p  at  y  =  rr,f  =  2.5. 


4.5  4.6  4.7  4.8  4.9  5  5.1  5.2  5.3  5.4  5.5 

Figure  5(1.  Coinpari.son  of  the  imnierieal  results  of  e  in  the  cusp  region  for  two 
different  numerical  resolutions:  512^  and  1024^. 
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5P10,256' 


.S’P10,512^ 


5P10. 1024^ 
ENOi.^Vi^ 


Figure  6a.  Time  history  of  the  maximum  density  from  different  computations. 
From  top  to  bottom:  SPlO  on  the  256^  grid,  SPlO  on  the  512^  grid,  SPlO  on  the  1024^ 
grid  and  EN03  on  the  512^  grid. 


EiV03,512* 


5P10, 1024* 


5P10,512* 


5P10, 256* 


Figure  6b.  Time  history  of  the  minimum  density  from  different  computations. 
From  bottom  to  top:  SPlO  on  the  256*  grid,  SPlO  on  the  512*  grid,  SPlO  on  the  1024* 
grid  and  EN03  on  the  512*  grid. 


Figure  9b.  Density  contour  at  *  =  '  ’for  the  initial  data  (4.4). 
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Time  history  of  u,x  at  the  cap 


Figure  9c.  Time  history  of  the  strain  rate  at  the  cap. 


Figure  10a.  Density  contour  at  t  =  3.5  for  another  set  of  initial  data. 
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.2 


Figure  10c.  Time  history  of  the  strain  rate  at  the  cap. 
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Figure  11.  Contour  plot  of  density  in  smother  set  of  calculation 
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